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Review  of  Some  Approximation  Operators  for 
the  Numerical  Analysis  of  Spectral  Methods 


Yvon  Maday 


Abstract.  This  paper  reviews  some  operators  that  are  used  in  the  nu- 
merical analysis  of  spectral  and  spectral  element  methods.  We  motivate 
the  introduction  of  these  different  operators  and  sketch  their  approxima- 
tion properties.  Finally,  we  apply  them  to  derive  optimal  error  estimates 
for  spectral  type  approximations  of  the  solution  of  elliptic  partial  differen- 
tial equations. 


§1.  Introduction 

Spectral  type  methods  are  high  order  discretizations  that  allow  to  compute 
approximate  solutions  of  partial  differential  equations.  The  recent  version  of 
spectral  approximations  is  based  on  the  Galerkin  approach  where  the  varia- 
tional statement  (equivalent  to  the  strong  formulation  of  the  PDE)  is  set  on 
discrete  spaces  of  test  and  trial  functions.  For  instance,  let  us  consider  the 
problem:  find  u G X such  that 

a(u,v)  = (f,v),  VuGX,  (1) 

where  X is  some  Hilbert  space,  and  a is  a continuous  bilinear  form  over  X. 
The  general  Galerkin  approximation  of  this  problem  first  requires  the  choice 
of  a family  of  discrete  spaces  X^  C X,  where  N is  a parameter  that  tends  to 
infinity  and  is  related  to  the  dimension  of  the  discrete  space  Xn-  The  discrete 
problem  is  then  stated  as  follows:  find  un  G Xjv  such  that 

u(un,vn)  — (/,ujv),  Vuat  G -Xjv-  (2) 

The  basic  general  hypothesis  that  makes  problem  (1)  well-posed  is  that  a is 
continuous  and  a-elliptic  over  X (i.e.  3 a > 0 such  that  a(u,  u)  > a||u||3f  for 
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all  u £ X.  These  properties  remain  true  over  each  Xn  (since  X,v  C X)\  thus 
(2)  is  also  well-posed  for  each  N.  In  addition,  the  solution  rty  satisfies 

||«  - WAr||x  < c inf  ||'U-vw||x-  (3) 

The  constant  c that  appears  in  (3)  is  the  quotient  of  the  continuity  constant 
of  a with  the  ellipticipty  constant  a and  is  thus  independent  of  Xpf. 

Going  back  to  spectral  methods,  the  definition  of  Xn  involves  polyno- 
mials, and  in  the  most  simple  cases  (we  shall  see  more  general  examples  in 
Section  5)  we  have  Xn  = X fl  P;y , where  Py  represents  the  set  of  all  poly- 
nomials of  (partial)  degree  less  than  or  equal  to  N . Here  N is  the  parameter 
responsible  for  the  convergence  of  the  method.  Due  to  (3),  one  ingredient  in 
the  numerical  analysis  of  the  spectral  method  is  the  approximation  properties 
of  the  space  of  polynomials  for  given  functions.  The  classical  analysis  of  the 
approximation  properties  of  polynomials  is  done  in  terms  of  L°°-  norms.  This 
is  not  completely  appropriate  for  our  purpose  since  most  often  X is  a Hilbert 
space  (generally  L2  or  H 1 spaces),  and  the  approximation  properties  have  to 
be  measured  with  these  norms.  If  a rate  of  convergence  (with  respect  to  N) 
on  the  best  fit  inf„N6xN  ||M  — uw||x  is  sought  after,  some  regularity  has  to  be 
assumed  over  u.  In  Section  2,  we  give  a survey  of  these  best  approximation 
results  depending  on  the  regularity  of  the  function  we  want  to  approximate. 
We  first  analyze  the  L2  best  fit  and  then  the  if 1 -best  fit.  The  main  ingredient 
in  this  analysis  relies  on  the  Legendre  basis  that  is  composed  of  the  orthogo- 
nal polynomials  for  the  standard  Lebesgue  mesure  over  the  interval  (—1,  +1). 
These  polynomials,  denoted  as  (Ln)n,  are  defined  by:  degree(Ln)  = n, 

Ml)  = 1,  (4) 

= (5) 

They  satisfy  some  standard  properties  (actually  valid  for  most  families  of 
orthogonal  polynomials) 

A(Ln)  = -±((l-?)^=n(n  + l)Ln,  (6) 

that  one  can  translate  by  saying  that  the  Legendre  polynomials  are  the  eigen- 
vectors of  the  (Sturm-Liouville)  operator  A.  Since  this  is  a possible  basis  set 
for  the  implementation  of  problem  (2),  this  gives  the  name  of  spectral  to  the 
methods  we  shall  consider  hereafter,  and  that  have  been  first  analyzed  in  [10]. 
We  refer  also  to  [6]  and  [3]  for  more  recent  surveys  on  the  numerical  analysis 
of  these  methods.  In  Section  3,  we  introduce  the  notion  of  numerical  integra- 
tion and  the  interpolation  operator,  two  notions  that  are  naturally  quite  close 
and  that  allow  to  transform  the  “theoretical”  approximation  method  into  a 
“applicable”  one.  In  Section  4,  motivated  by  the  analysis  of  the  Stokes  prob- 
lem, we  introduce  a new  operator,  that,  in  opposition  to  the  previous  ones,  is 
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uniformly  stable  (in  N)  both  in  the  i2-norm  and  the  if1 -norm  and  possess 
optimal  approximation  properties.  It  has  to  be  said,  beforehand,  that  in  the 
precise  analysis  of  these  spectral  (or  polynomial)  approximation,  the  Bern- 
stein inequality  runs  counter  to  most  standard  tools  that  generally  allow  for 
deriving  approximation  results  for  a new  operator  from  an  already  analyzed 
one.  This  Bernstein  inequality  tells  about  the  equivalence  of  norms  on  the 
finite  dimensional  linear  space  of  polynomials.  It  is  well  known  that,  for  any 
function  in  H1,  the  L 2 norm  is  smaller  than  the  If1 -norm;  of  course  this  is 
true  in  particular  on  polynomials: 

W>JV  € Pjv,  ||$n||l*  < 

Since  all  norms  are  equivalent  on  Pat,  there  exists  a constant  (obviously  de- 
pending on  N)  such  that 

V</>at  £ Pat,  W^nWh1  < c(A0||</>at||l2- 

The  behaviour  of  this  constant  is  made  precise  by  the  Bernstein  inequality 

V0jv  € Pat,  || (PnWh1  < cN2\\(I>n\\L2, 

where  c no  longer  depends  on  N.  This  estimate  is  optimal  (in  the  sense  that 
there  exists  a sequence  of  polynomials  such  that  the  ratio  of  the  II 1 -norm  over 
the  L2-norm  scales  like  0(N2)),  but  is  bad  as  regards  the  ratio  of  convergence 
rate  between  the  H1- best  fit  and  the  L2-best  fit  that  scales  like  C*(iV_1),  as 
we  shall  see  below. 

In  the  first  three  sections,  the  domains  where  the  functions  live  will  be 
very  simple,  actually  too  simple  to  tackle  real  life  problems;  indeed  these  are 
bricks  equal  to  (— l,l)d  where  d = 1,2  or  3.  The  generalization  of  spectral 
methods  to  more  complex  geometries  is  done  by  combining  two  key  ingredi- 
ents: the  mapping  of  bricks  onto  curved  bricks  through  regular  mappings,  and 
domain  decomposition.  We  give  some  hints  about  this  generalization  in  §5. 

§2.  Hilbert  Type  Projection  Operators 

Let  us  start  with  the  one-dimensional  case.  In  L2(— 1, 1),  we  consider  the  set 
Pat(— 1, 1)  of  all  polynomials  of  degree  < N,  From  the  Weierstrass  density 
theorem,  we  know  that  any  element  <j>  in  L2(— 1, 1)  can  be  written  as 

OO 

*(C)  = £?"£»(0,  (7) 

n=0 


where  the  convergence  of  the  series  holds  in  L2.  The  coefficients  <f>n  can  be 
derived  from  (j>  thanks  to  the  orthogonality  of  the  Legendre  basis  as  follows: 


4>n  = 


2 n -j- 1 


j\(()Ln«)d(. 


2 
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Next,  from  (6)  we  derive  that 


f ♦k)41 

J- 1 n(n 


ALn{ 0 

+ 1) 


dC, 


noticing  that  „4  is  symmetric,  and  assuming  <j>  regular  enough,  we  derive  that 
2n  + 1 


^n  = 


• t Amo-f^dc 

7-i  w(n  + 1) 


If  we  iterate  this  argument  p times,  we  obtain 


2n  + 1 


['  Ap(<t>){C) d(, 

7_ i ^ASV(n  + l)P  s’ 


so  that,  the  following  simple  relation  holds  between  the  Legendre  coefficients 
of  4>  and  of  Ap{<j>)'. 

= , 1 M4>)n ■ 

nP(n  + 1)p 

Next,  let  7 tn  denote  the  L2  ( — 1 , l)-projection  over  P^(— 1, 1).  Going  back  to 
(6),  we  deduce  from  (7)  and  (5)  that 


N 

*„(<£)  = 0,  (8) 

n=0 

so  that 

OO  OO  1 

E *"MC)  = E nPfn+np-4^"^ 

n=7V+l  n=JV+l  ' ^ ’ 

and,  by  Parseval 

oo  1 „ 

n=lV+l  ' ' 

5 ivi4*  £ 

n=AT+l 

1 °°  o 1 

s i^r  - [Fni^w)in.,-u,' 

n=0 

We  have  thus  proven  that,  for  any  (j>  in  the  domain  V[AP]  of  _4P, 

\\<t>  - < c(p)N-2P||^P(^)||i2(_1,1). 

It  is  easy  to  check  that  H2p(— 1, 1)  C V\AV]\  hence  the  following  theorem  (due 
to  Canuto  and  Quarteroni  [7]),  proven  here  for  even  values  of  r,  holds  for  any 
r thanks  to  an  argument  of  interpolation  between  Sobolev  spaces: 
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Theorem  1.  For  any  real  number  r > 0,  there  exists  a constant  c > 0, 
depending  only  on  r such  that,  for  any  function  cf>  € HT(— 1, 1), 

||<£  - 7fiv(0)||z,2(— i,i)  ^ ciV-r||0||JJr(_i,1).  (9) 

Let  us  denote  now  by  II  jv  the  L2((— 1,  l)d)  orthogonal  projection  operator 
over  the  set  Pat((— l,l)d)  of  all  polynomials  of  degree  < N with  respect 
to  each  variable.  By  Fubini’s  theorem,  IIjv  = ttn  ® ttn,  in  2D  and  fJw  = 
7 tn  ® 07 tjv  <8>  o7T;v,  in  3D.  By  tensorizing  (9)  we  derive 

Theorem  2.  For  any  real  number  r > 0,  there  exists  a constant  c > 0, 
depending  only  on  r such  that,  for  any  function  0 € ffr((—l,  l)d), 

||0  - niV(0)IU2((-1,1)a)  < cN-rM\\Hr{(_hl)d). 

We  are  now  in  a position  to  tackle  the  approximation  in  the  H 1 norms. 
First,  we  consider  a function  0 € Hq (— 1,  l)fl/fr(— 1, 1),  with  r > 1.  It  is  quite 
immediate  to  check  that  the  polynomial  0jv(C)  = fl ] 11 N belongs 
to  Pat(— 1,1),  vanishes  at  £ = —1,  and  satisfies 

0n(1)  = J = J 7rAr-i^(£)-M0rf£ 

= 0(1)  - 0(-l)  = 0, 

and  hence  is  an  element  of  Pjv(— 1, 1)  (1  1).  Finally  it  is  a good  ap- 

proximation of  0,  since  from  Poincarre’s  inequality  and  Theorem  1, 

110  ~ 0iv||Hi(_l,l)  <C\\~-  lil.2(  — 1,1) 

< cAr1-r||^||Kr-l(_lil)  < dV1-10||^(-1,l). 

Let  us  introduce  now  the  orthogonal  projection  operator  7rjj°  from  Hq(— 1, 1) 
onto  Pat(— 1, 1)  fl  Hq(— 1, 1),  we  can  state  the  following  result  (due  to  Maday 
and  Quarteroni  [15]): 

Theorem  3.  For  any  real  number  r > 1 and  any  real  number  0 < s < 1, 
there  exists  a constant  c > 0,  depending  only  on  r and  s such  that  for  any 
function  0 £ 1)  fl  Hr(~  1, 1), 

110  - 7Lv°(0)lltfq-i,i)  ^ cNs~r  j|0||rrT'(— i,i)*  (10) 

Proof:  The  theorem  has  been  obtained  for  s = 1.  For  s = 0 it  is  obtained 
through  a standard  Aubin-Nitsche  duality  argument,  and  then  for  any  s by 
interpolation  between  Sobolev  spaces.  □ 
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Remark.  At  this  point  it  has  to  be  said  that  the  L 2 projection  operator  ttn 
does  not  have  optimal  approximation  properties  in  the  //'--norm,  the  only 
(non-improvable)  property  that  can  be  obtained  is 

\\<t>  - 7rjv(<wil/p(-i,i)  < cN^~r\\(f>\\Hr{_hl). 

We  refer  to  [3]  for  details  and  counter-examples. 

Remark.  It  may  also  be  interesting  to  note  that,  despite  their  definition,  the 
previous  operators  have  stability  properties  in  various  norms.  First  for  the 
L2-operator,  we  have 


||7TjV0||//i(-l,l)  < ciVa||<£||tfi(_lil), 

which  is  related  to  what  we  have  indicated  in  the  previous  remark,  but  also 

lkwVlU2(-i,i)  < 

which  is  rather  suprising  since,  from  this  (non-uniform)  stability  property,  the 
Hq -projection  operator  can  be  extended  to  (irregular)  functions  of  L2!! 

Again  by  tensorization  of  the  results  of  the  one  dimensional  Theorem  3, 
we  exhibit  a polynomial  that  approximates  regular  functions  in  JTq((-1,  l)d) 
well,  from  which  we  derive  approximation  properties  on  the  multidimensional 
projection  operator  Ilj^0  from  R^((-l,  l)d)  over  Pw((-1,  l)d)  n Rg((-1,  l)d) 
both  in  R1-norm  and  in  L2-norm  (derived  by  duality): 

Theorem  4.  For  any  real  number  r > 1 and  any  real  number  0 < s < 1, 
there  exists  a constant  c > 0,  depending  only  on  r and  s such  that,  for  any 
function  <f>  € (( — 1,  l)d)  n Hr((- 1,  l)d), 


||0  - n^°(</i)||//S((_lil)d)  < cAs  T'||0||//r((— (if) 

These  results  can  be  completed  in  order  to  derive  a whole  scale  of  ap- 
proximation projectors  in  higher  order  norms.  These  are  required,  e.g.  for 
the  analysis  of  the  approximation  of  fourth-order  problems.  The  general  re- 
sult, concerning  the  orthogonal  projection  operator  II from  // p ( ( — 1 , l)d  n 
Hq  ((— 1,  l)d  onto  Pyy((— 1,  l)d)  fl  Hq  ((—1,  l)d)  is  given  in  the  following  theo- 
rem (due  to  Maday  [11]  in  ID,  see  also  [3]  for  the  extension  to  2 and  3D): 

Theorem  5.  For  any  real  number  0 < a < p and  any  0 < s < p < r,  there 
exists  a constant  c > 0,  depending  only  on  r,  s,  p,  a such  that,  for  any  function 

0gRoct((- i,i)d)nF((-i,i)d), 


H°(( 


-l.i)**)  < CNS 
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Remark.  A final  remark  on  these  operators  is  that  improved-approximation 
results  in  negative  norms  are  also  true,  and  can  be  obtained  in  a classical  way, 
by  further  refering  to  the  Aubin-Nitsche  duality  argument.  Hence,  Theorem 
5 is  also  valid  for  negative  values  of  s. 

These  results  allow  us  to  prove  that  the  approximation  of  most  elliptic 
variational  problems  by  spectral  methods  is  optimal.  As  an  example,  let  us 
consider  the  (non-constant)  Laplace  problem  on  a cube  ft  = ( — 1,  l)3:  given 
a 3 x 3 matrix,  symmetric  and  uniformly  positive  definite,  we  consider  the 
problem  of  finding  u G Hq  (ft)  such  that 

— div[Agrad]u  = f.  (12) 

The  approximation  then  consists  in  finding  an  element  un  in  = Pjv(ft)  fl 
Hq  (ft)  such  that 


/ AVwjvVu n = [ fvN i e Xn.  (13) 

Jci  Jn 


Assuming  that  u € i/r(ft),  we  deduce  from  (3)  and  (11)  that 


||m  - < cN1  r|M|ff’-(n). 


As  hinted  in  the  introduction,  this  problem  is  numerically  intractable; 
indeed  the  implementation  of  (13)  requires  the  computation  of  the  two  in- 
tegrals appearing  on  the  left-  and  the  right-hand  sides  of  this  equation.  The 
exact  computation  is  most  often  impossible,  and  certainly  numerically  not  fast 
enough.  The  use  of  numerical  integration  rules  is  the  cure  to  this  problem, 
but  in  order  to  combine  efficiency  and  precision,  following  Gottlieb  [9]  and 
Mercier  [17],  we  refer  to  the  use  of  Gauss  type  quadrature  rule.  Indeed,  they 
are  well  known  to  be  well  suited  for  the  integration  of  polynomials. 

§3.  Interpolation  Operators 

Between  the  different  numerical  quadrature  rules  over  (—1,1),  well  suited  for 
polynomial  integration,  we  shall  quote  here  the  Legendre- Gauss  and  Legendre- 
Gauss-Lobatto  ones.  We  refer  to  [2]  for  more  details.  For  the  sake  of  com- 
pleteness, we  recall  the  definition  of  these  formulae: 

Theorem  6.  (Gauss  formula)  For  any  real  number  n,  there  exists  a unique 
set  of  points  — 1 < < Q < ■ ■ • < < 1>  and  a unique  set  of  positive 

weights  u>”  such  that  for  any  polynomial  cj>  € P2n_i(— 1,1),  the  following 
equality  holds: 

f md<  = i>« ?k- 

-'-i  i=i 
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Theorem  7.  (Gauss-Lobatto  formula)  For  any  real  number  n,  there  exists 
a unique  set  of  points  — 1 = < £?  < • • • < £"  = 1,  and  a unique  set 

of  positive  weights  p"  such  that  for  any  polynomial  (j>  € P2n-i(— 1,1),  the 
following  equality  holds: 

f1  mdc  = j2m)p?- 

»= 0 

Prom  now  on,  we  shall  assume  that  the  degree  of  the  polynomials  for  the 
approximation  is  fixed  to  be  N,  and  we  shall  use  N + 1 points  either  of  Gauss 
or  Gauss-Lobatto  type.  For  the  sake  of  simplicity,  these  points  will  be  denoted 
with  no  superscript,  i.e.  in  all  of  what  follows,  we  set  Q = C,N+1  and  £,  = . 

We  recall  that  these  points  are  the  roots  (resp.  the  extrema)  of  the  Legendre 
polynomials;  more  precisely,  we  have 

Vi,  Ln+ i(Ct)  = o,  (resp.  (1  - g)L'N+1(Zi)  = 0). 

After  tensorization,  these  one  dimensional  quadrature  rules  easily  provide 
quadrature  rules  on  the  square  and  on  the  cube  defined  as  follows  (e.g.  in 
2D  for  the  Gauss  Lobatto  formula): 

N N 

GL  i=0j=0 

The  problem  that  is  actually  implemented  is  then  the  following:  find  an  ele- 
ment «jv  in  Xn  such  that 

y^AVu^yvN  = y ^fvN,  Vv/v  e Xn.  (14) 

GL  GL 

Even  in  the  case  where  A is  constant,  at  least  in  more  than  one  dimension, 
the  left-hand  side  is  not  exactly  computed.  The  problem  is  no  longer  of  the 
form  (1),  and  the  abstract  theory  has  to  be  generalized  in  order  to  handle  this 
problem  as  well. 

Here  is  not  the  place  to  detail  this  generalization  (see  [3],  where  the 
complete  analysis  is  performed)  but  it  is  natural  that  the  a-ellipticity  of  the 
bilinear  form  on  the  left-hand  side  of  (14)  is  again  one  of  the  key  ingredients 
and  has  to  be  satisfied.  This  follows  from  the  property,  proven  in  [7] 

W>jv  G Pjv(-1,  1),  X>*  > f <t>2N(Od(. 

GL  J~1 

From  this  property  it  can  be  easily  derived  that  the  solution  u ^ to  (14)  exists 
and  is  unique. 

The  approximation  properties  of  the  polynomial  interpolation  operator 
over  the  Gauss-Lobatto  nodes  is  of  great  importance  in  the  error  bounds.  Let 
jjv  denote  this  operator  in  one  dimension: 

V<£eC°([-l,l]),  «;v(0)  £ Pjv(— 1, 1)  and  Vi,0  < i < N,  = <K&) 

and  let  us  tensorize  it  in  order  to  get  a two  (resp.  a three)  dimensional  operator 
IN  = in  ® in  (resp.  1 n = in  ® in  ® in).  The  properties  of  this  operator  have 
been  established  in  [12]  and  [2],  and  read  as  follows: 
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Theorem  8.  For  any  real  numbers  s and  r satisfying  r > (d  + s)/ 2 and 
0 < s < 1,  there  exists  a positive  constant  c depending  only  on  r such  that 
for  any  0 in  Hr((—1,  l)d)  the  following  estimate  holds 

110  -2jv($)||ff»((-i,i)<«)  < cArs_r||0||/f((_iii)d).  (15) 

It  has  to  be  noticed  that  this  operator  requires  more  regularity  than  the 
L2  projection  operator,  but  it  is  optimal  both  in  the  L 2 and  the  H1  norms. 
It  has  also  to  be  recalled  that  in  the  classical  approximation  properties  in  the 
L°°  norm,  the  Lebesgue  constant  appears  as  a pollution  of  the  approximation 
properties  of  the  interpolation  operator  as  regards  the  optimality  provided 
by  the  corresponding  best  fit.  This  is  not  the  case  in  the  L2-norm.  In  this 
direction  what  we  have  more  precisely  is  that,  for  any  function  0 in  Hq(— 1, 1), 

ll*w0IUa(-i,i)  ^ c(ll0IUa(-i,i)  + ]v  II^Hi2(— i-i))> 

and  for  any  function  0 in  FF1  ( — 1, 1), 

||*JV0||ifl(-l,l)  < C||0||ffl(-l,l)- 

Another  nice  property  of  this  operator,  that  has  some  importance  for  nonlinear 
PDE’s,  is  the  following  result:  for  any  polynomial  0m  € Pm(-1,  1), 

M 

||»N0m||l2(-1,1)  < c(  1 + -jy)||0Ar|Ua(-l,l)- 

Here  no  duality  argument  allows  us  to  derive  from  the  previous  theorem 
improved  approximation  properties  in  negative  norms.  It  is  an  open  problem 
to  derive  such  results. 

The  numerical  analysis  of  problem  (13)  then  continues  by  noticing  that 

Yjvn  = YiN(f)vN, 

GL  GL 

which  is  one  of  the  ingredients  that  allows  to  prove  (see  [2]): 

Theorem  9.  Assume  that  the  solution  u of  (12)  belongs  to  Hr(fT),  that  the 
coefficients  in  A are  very  regular,  and  that  the  data  f belongs  to  Then 

the  solution  un  to  (13)  satisfies 


||«  - «/v||ifi(n)  < c(-^1  7lM|j/’-(n)  + N Pll/Ilfl>(n))- 

The  case  where  A is  not  so  regular  can  be  handled  with  the  same  type 
of  arguments,  but  more  technical  tools  are  involved;  we  refer  to  [16]  for  more 
details.  It  is  interesting  also  to  note  at  this  level  that,  taking  into  account  non- 
homogeneous  Dirichlet  boundary  condition  is  very  simple  thanks  to  the  nice 
properties  of  the  interpolation  operator.  Indeed,  assume  that  the  solution  to 
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our  problem  (12)  has  to  satisfy  (instead  of  zero  Dirichlet  boundary  conditions) 
the  following  condition:  u\qq  = g where  g is  a given  function  on  the  boundary 
of  ft.  Then,  naturally,  for  the  approximation,  we  look  for  a polynomial  un  in 
Piv(ft)  such  that  (12)  holds,  and  in  addition 

= T-nQ, 

where  In  is  the  operator  of  interpolation  defined  edge  by  edge  (respectively 
face  by  face)  from  iyv  (resp.  from  In)-  Since  the  interpolation  operator  is 
optimal  both  in  L2  and  in  H1 , it  results  by  an  argument  of  interpolation 
between  Sobolev  spaces  that  it  is  also  optimal  with  respect  to  the  Hl/2(dfl)~ 
norm.  This  fractional  order  norm  is  the  natural  one  for  the  treatment  of  the 
boundary  terms.  It  has  also  to  be  stressed  that  neither  the  L2  projection 
operator  nor  the  Lf1 -projection  operator  allow  such  an  optimality  nor  such 
ease  in  the  implementation. 

Next,  associated  with  the  Gauss  quadrature  formula,  we  can  also  define  an 
interpolation  operator,  denoted  as  j,x  and  defined  as  follows:  V</>  6 C°([—  1, 1]), 

jN{4>)  € Pat+i(-1,  1)  and  Vi,  1 < i < N + 1,  jjv^XCi)  = 

The  L2{— 1,  ^-approximation  properties  of  this  second  interpolation  operator 
are  also  optimal.  Unfortunately,  in  the  H 1 -norm  it  is  not  optimal;  for  instance 
it  is  readily  checked  that  ~ Ln-i)  — Ln-i-  Recalling  that 

Vn,  £„+!-!„_!  = — 2-^\-(l-C2)Lln, 

nyn  + 1) 

it  is  then  easily  proven  that  ||Ljv+i  — Z/^v— 1||//1(— i,i)  scales  like  C>(Arl|/2)  while 
||Ljv-i||i/i(-i,i)  scales  like  0{N)\  jn  is  thus  not  stable  in  the  H 1 norm. 

For  similar  reasons,  the  interpolation  operator  In  on  the  Gauss-Lobatto 
nodes  does  not  have  optimal  approximation  properties  in  the  H2(—  1,  l)-norm. 
In  order  to  achieve  such  a property,  we  have  to  refer  to  generalized  Gauss- 
Lobatto  rules  as  is  done  e.g.in  [1], 

§4.  An  “Ideal”  Operator 

At  this  stage  there  is  no  operator  from  L2(— 1,1)  onto  the  set  of  polynomials 
that  has  optimal  approximation  properties  and  is  stable  both  in  the  L2  and 
the  H 1 norm.  Such  an  operator  is  useful,  as  will  be  explained  below,  in  the 
analysis  of  the  Stokes  problem.  In  order  to  define  this  “ideal”  operator,  we 
fix  a positive  real  number  A and  a cut-off  function  x of  class  C1  on  1R+  such 
that  x is  equal  to  1 on  [0, 1 — A],  decreases  from  1 to  0 on  [1  — A,  1]  and 
vanishes  on  [l,oo].  Next,  with  each  positive  integer  N , we  associate  as  in  [18] 
an  operator  with  values  in  P/v(— 1, 1)  fl  Hq(—1,  1)  as  follows:  since  each 
function  <f>  in  Hq(— 1, 1)  can  be  written  as  </>  = 4>n{Ln+\  — Ln_ i),  we  set 

x(jr)4>n {Ln+ 1 — Ln- 1).  Note  that  the  sum  above  is  finite  since 
X has  a bounded  support.  It  is  proven  in  [4]  that  this  operator  is  stable  both 
in  the  Hq  and  the  L2  norms: 
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Theorem  10.  There  exists  a constant  c,  independent  of  N such  that  for  any 
function  </>  € Hq(-  1, 1), 


ll7r&0lUJ(-l,l)  ^ cll0IUa(-l,l)>  ll7rAT^II w1  ( — 1,1)  ^ II^IIhH-1,1)-  (16) 

It  is  an  easy  matter  to  verify  that  the  operator  leaves  invariant  all 
polynomials  of  Pajv(~  1, 1)  fl  Hq(— 1,  !)•  The  previous  stability  and  the  best 
fit  estimates  (9),  (10)  imply 

Theorem  11.  For  any  positive  real  number  r and  any  real  number  0 < s < r, 
there  exists  a constant  c > 0,  depending  only  on  r and  s such  that,  for  any 
function  <j>  6 Hq(—1,  1)  if  r < 1 and  any  function  (f>  € 1)  fl  Hr(—1, 1) 

if  r > 1, 

\\<t>  — 7rTV'^l|ifs(— 1,1)  5:  cNs~r  ||<^j|ff(_l,i). 

As  an  application  of  the  previous  result,  we  can  consider  the  problem 
of  finding  compatible  spaces  for  the  approximation  of  the  Stokes  equation. 
Under  variational  formulation,  this  problem  consists  in  finding  a pair  (u,  p)  in 
(.ffo(fl))d  x Lo(fi)  of  velocity  and  pressure  such  that 


/ 
J u 


VuVv  - 


fv,  vt ,e(Hj(n))d, 
Vg  e L§(ft), 


(17) 

(18) 


where  Lq(U)  *s  the  set  of  L2  functions  with  zero  average.  It  is  well  understood 
now  that  the  spectral  approximation  of  the  Stokes  problem  based  on  polyno- 
mials of  the  same  degree  leads  to  instabilities.  This  is  due  to  the  fact  that 
the  pressure  space  is  too  rich  in  comparison  to  the  velocity  space.  Indeed, 
there  exist  polynomials  gjv  in  Pjv((— 1,  l)d)  such  that  J^qN  divv/v  = 0 for  all 
vN  in  (Pw((-l,l)d)n.ff01((-l,l)d))d  (e.g.  qN(x,y,z ) = LN(x)LN(y)LN(z)). 
Of  course  such  polynomials  (called  spurious  modes)  prevent  the  discrete  prob- 
lem from  being  well-posed  since  it  prevents  the  definition  of  a unique  pres- 
sure. The  cure  is  well  known,  and  consists  in  depleting  the  pressure  space 
for  a given  velocity  space.  In  [14]  the  pair  (Pat((— 1,  l)d)  fl  Hq((— 1,  l)d))d  x 
PJV_a((-l,  l)d)nLg((-l,  l)d)  has  been  proposed,  and  gets  rid  of  the  spurious 
modes.  It  is  known  as  the  Pat  x Pjv_2-method.  Actually,  what  is  looked  for 
is  a pair  XN  x M/v  approximating  1,  l)d))d  x Lq((— 1,  l)d)  well  and 

such  that  not  only  V<pv  € Mjv,3u/v,  fnpdivvj^  / 0,  but  more  precisely,  in 
order  to  get  a stable  method,  we  require  that 

Vg/v  € Mn,3vn,  / pdiv«AT  > /?||vJv|lfl-i((_i,iW)||gAr||i2((_i>1)d), 

Jti 


where  ,6  is  known  as  the  constant  of  the  inf-sup  condition.  The  behaviour 
of  (3  for  the  P.v  x Pjv-2~method  scales  as  0(N~~ ~)  (see  [2]),  and  it  has 
been  a long  standing  question  whether  there  is  a uniformly  stable  spectral 
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approximation  of  the  Stokes  problem.  It  has  to  be  said  that  the  nonuniform 
behaviour  of  the  inf-sup  constant  pollutes  the  accuracy  of  the  pressure,  but 
also  pollutes  the  convergence  properties  of  some  classical  solvers  for  the  Stokes 
problem  (see  [13]).  The  “ideal”  operator  introduced  above  allows  us  to  prove 
that  a compatible  choice  is  the  Pa?  x P^A'-method  that  proposes,  for  the 
same  choice  of  velocity  space,  Paa?((~  1>  l)d)nLo((— 1,  l)d)  to  be  the  pressure 
space.  The  following  result  is  due  to  Bernardi  and  Maday  [4] : 

Theorem  12.  For  any  real  number  A,  0 < A < 1,  there  exists  a positive 
constant  (3  independent  of  N such  that,  for  any  integer  N > 2/(1  — A)  and 
any  qN  € Paa?((-1,  l)d)  fl  Lg((- 1,  l)d), 


sup  Jn  pdivUN 

tfwe(Pw((-i,i)d)ntf>((-i, !)<<))“  l|vAr||//i((_i,i)<*) 


> /3||<?A?||z,2((-l, !)■<)• 


Proof:  Let  <7a?  be  any  polynomial  in  Paat((— 1,  l)d)  fl  Lq((— 1,  l)d).  It  is  a 
standard  matter  (see  e.g.  Corollary  2.4  in  [8])  that,  to  q a?,  can  be  associated 
a (continuous)  element  v in  [/fg((— 1,  l)d)]d  such  that 


divu  = qN  and  |M|tfi((_i,i)<<)  < c||gN||x,2((-i,i)d)- 

The  problem  is  that  v is  not  a polynomial.  We  define  «a?  = ® *n  2D 

and  v a?  = ^ at  ® ^ a?  ® ttJ/v  'n  3D  ^or  wh'ch  we  derive  thanks  to  (13)  that 

ll«Ar||j/i((-i,i)d)  < cIIq,a? Ilr,2((— i,i)d)- 

Due  to  the  fact  that  7r^  leaves  invariant  all  polynomials  of  Paa'(— L 1)  H 
Hq(— 1, 1),  we  deduce  that  fn  qwdiv(vw  — v)  = 0,  and  thus 

/ qNdivvN  = / qN  divv  = / q2N, 

Jn  Jn  Jn 

which  concludes  the  proof  with  /?=!.□ 


§5.  Extension  to  Domain  Decompositions 

In  the  spectral  method  history,  the  need  to  tackle  more  general  domains  was 
recognized  early.  In  this  direction,  Patera  has  proposed  in  [19]  the  spectral 
element  method  that  combines  the  accuracy  of  the  spectral  method  with  the 
flexibility  of  the  domain  decomposition  methods.  The  idea  is  to  introduce  a 
partition  of  the  domain  ft  as  a union  of  nonoverlapping  subdomains: 

ft  = u£=1ft\  ft*  n ftf  = 0. 

In  addition,  we  assume  that  each  subdomain  ft*  is  associated  with  a regular 
one-to-one  mapping  Tk  that  maps  the  brick  (— l,l)d  onto  ft*  and,  for  the 
time  being  at  least,  we  make  the  following  assumptions: 
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Assumption  1. 

' an  entire  common  face  (in  3D),  or 
an  entire  common  edge,  or 

a common  vertex,  or 

u 

Assumption  2.  The  two  parametrizations  of  the  previous  intersection 
^ £ 

O Dfi  , resulting  from  Tk  and  coincide. 

This  allows  us  to  define  the  discrete  space 

XN  = {vN  6 vN[n,  o Tk  € Pw((-1,  l)d)} 

and  the  discrete  associated  problem  (2)  (or  its  implementable  version  involving 
the  Gauss-Lobatto  quadrature  rule  over  each  ftfc  as  in  (13)). 

The  main  ingredient  that  allows  us  to  prove  that  the  previous  scheme  is 
again  optimal  lies  in  the  definition  of  a element  in  Xu  that  approximates  well 
a given  regular  function  u.  This  is  done  easily  by  considering  the  element  ujv, 
defined  locally  over  each  subdomain  as  oTk  = Jjv[w|n*  It  results 

from  Assumptions  1 and  2 that  vn  is  actually  continuous  and  vanishes  over 
dft.  From  (15)  it  is  an  optimal  approximation  of  u in  the  sense  that 

11“  - vjvIlH^n)  < ciV1-r||u||tf.-(n).  (19) 

The  best  fit  in  H1  (12)  is  certainly  as  good  as  the  proposed  vn,  and  the  spectral 
element  method  can  be  proven  to  be  an  optimal  approximation.  We  have  only 
sketched  the  numerical  analysis  of  this  approximation,  since  the  main  purpose 
of  this  paper  is  to  discuss  projection  operators.  It  is  fundamental  to  have  used 
here  the  interpolation  operator  to  construct  v/v,  since  it  provides  a globally 
continuous  function.  As  an  example,  the  use  of  the  //^-projection  operator 
would  not  have  given  rise  to  a continuous  function  since,  for  a given  function 
<j>  over  the  brick  (—1,  l)d,  the  value  of  njy°(0)  over  any  face  depends  not  only 
on  the  value  of  4>  on  the  given  face,  but  depends  on  cj>  inside  the  whole  domain. 

We  want  to  end  this  section  by  giving  some  hints  on  the  “mortar  spectral 
element  method”  due  to  Bernardi,  Maday  and  Patera,  that  allows  to  relax 
assumptions  1 and  2 (and  even,  more  generally,  allows  to  combine  spectral 
methods  on  some  subdomains  with  different  finite  element  methods  on  others 
see  [5]).  Due  to  lack  of  space,  but  also  in  order  to  better  understand  the 
main  feature  of  the  projection  operators  that  is  at  the  basis  of  the  method, 
we  shall  consider  a simple  two  dimensional  domain  ft  = (—1,2)  x (-1,1) 
decomposed  into  3 subdomains  ft1  = (—1,1)  x (—1,1),  ft2  = (1,2)  x (—1,0) 
and  ft3  = (1,2)  x (0,1).  This  decomposition  violates  assumption  1 since  the 

l 2 

intersection  ft  fl  ft  is  not  a common  whole  edge.  We  want  nevertheless  to 
propose  a discrete  method  that  will  allow  to  provide  an  optimal  approximation 
of  the  solution  u of  (12)  (with  A =Id  for  the  sake  of  simplicity).  The  discrete 
space  X*N  that  we  propose  is  imbedded  in 

Yn  = {“n  € L2(ft),uyV|fj»;  € = 0,vn\q2  = “win3  over  fl  nft  } 


ft*  n ft*  = 
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but  it  is  readily  checked  that  imposing  continuity  at  the  level  of  the  inter- 
face x = 1 will  rigidify  the  approximation  and,  in  the  general  case,  will  spoil 
the  accuracy  of  the  method.  In  order  to  relax  this  continuity  condition  (re- 
mind that  it  is  inherited  from  the  requirement  that  X*N  C X),  we  resort  to 
nonconforming  approximations.  We  shall  replace  the  continuity  condition  by 
requiring  that,  over  the  interface  x = 1,  we  impose  for  each  element  in  Yn 

J [v~(l,y)  ~v+(l,y)]ipN{y)dy  = 0,  VipN  € Pat-2(-1,  1),  (20) 


where  v = U|ni  and 

v+  = / V\Q2  for  (x<y)  e 

( V|n3  for  ( x,y ) £ f l3. 

Since  v~(l,y)  has  to  vanish  for  y = ±1  (due  to  the  homogeneous  boundary 
conditions),  it  is  entirely  defined  by  the  N — 1 conditions  in  (20);  in  particular 
choosing  ip  in  Pjv(— 1,1)  would  be  much  too  stringent.  The  elements  of  Y)v 
that  satisfy  (20)  constitute  the  space  X*N  of  approximation.  The  method  is 
then:  find  UN  € such  that 


aN{vrN,vN)  = J2  [ Vu*nVvn  = v[  fvN,  VvN€X*N.  (21) 

k=iJ°k  k=iJnk 


Since  X*N  is  no  longer  a subspace  of  X,  the  ellipticity  of  the  bilinear  form  of 
this  problem  is  not  straightforward.  Nevertheless,  it  is  true  (and  here  it  is 
particularly  obvious  since  d£lk  fl  dft  / 0).  This  argument  allows  us  to  check 
that  there  exists  a unique  solution  u*N  to  (21).  In  order  to  derive  the  error 
bound  we  proceed  as  follows:  for  any  w/v  zx*N, 


a||«Ar  ~ wn II*  < aN(u*N  - wN,u*N  - wN) 
K „ K 


= yz  / vwnv(u*n~wn) 

k= iJnk  k= iJnt 

K f K [ 

= 'V'  / -Au(uf,~WN)~y2  VwnV(u*n-wn) 
k= iJnk  k^iJ*k 

= y2  VuV(u*N  - wN)  - Y]  VwjvVfu^-tojv) 

k= iJak  k= iJnk 

~ /.j  ^0^U*N  ~ WN ^ _ ^ ~ WN'>  + 


so  that,  from  (20)  we  derive  that  for  any  ip  £ P/v-2(— 1,1) 
K 


a\\uN  ~ WJV||* 


< V / X{u-wN)X(u*N  - 1VN) 

k=iJ*k 

- J ~ ^U*N  ~ WN')~  - (U*N  - Wjv^+5- 
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It  follows  from  the  previous  inequality  that 


a\\v,*N  - wjv||*  < c||u  - wjv||«  + sup 


ix=i[||  Z^]Kr  -vftdy 

IMI. 


(22) 


By  choosing  if>  equal  to  7r/v-2wr,  it  results  that 


sup  JU1&  *]K  v»]dy  <cN-r\\u\\Hr{n). 
vNex'  Pw||* 


It  remains  to  choose  a good  approximation  w jv  of  u in  X*N  to  take  into  account 
the  first  term  on  the  right  hand  side  of  (22).  This  is  done  by  noticing  that,  for 
any  <j>  £ Hq(—1,  1),  the  element  <f>N  of  Pjv(— 1, 1)  fl  Hq(— 1, 1)  that  satisfies 


[<t>N  - <f\i>N(y)dy 


= 0, 


VipN  S Pjv-a(-l,  l)i 


is  nothing  other  than  <j>u  = ir^0(4>).  Indeed,  we  remark  that,  for  any  xn  S 
Pat(-1,  1)  n #o(-l>  1),  then  Xn  € Pjv-2(-1,  1),  thus 

J [<I>n  - (f\xN(y)dy  = - / y>N  - <l>]'xN(y)dy- 

The  choice  of  a good  element  wn  is  done  as  follows.  We  first  set  wN\ nk  = 
In(un\cii>)  that  is  an  element  of  Ypf.  We  then  set  =WN\nk  for  A:  = 2, 3, 

and  build  wN\ni  by  adding  to  the  correction  7rjv°(®iv  — w] y){y)  ^ 
so  that  it  satisfies  (20).  Due  to  the  optimal  approximation  properties  of  the 
operator  both  in  the  L2  and  in  the  H1 -norms,  we  deduce  that  the  mortar 
spectral  approximation  (21)  is  optimal  in  the  sense  that  (19)  still  holds. 
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